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■ We present a semi-analytic model for the thermal and ionization history of the uni- 

verse at 1000 > z > 3. This model incorporates much of the essential physics included 

. in full-scale hydrodynamical simulations, such as (1) gravitational collapse and virial- 

ization; (2) star/quasar formation and subsequent ionizing radiation; (3) heating and 
cooling; (4) atomic and molecular physics of hydrogen; and (5) the feedback relation- 
ships between these processes. In addition, we model the process of reheating and 
reionization using two separate phases, self-consistently calculating the filling factor 
of each phase. Thus radiative transfer is treated more accurately than simulations to 
date have done: we allow to lowest order for both the inhomogeneity of the sources 
and the sinks of radiation. After calibrating and checking the results of this model 
against a hydrodynamical simulation, we apply our model to a variety of Gaussian and 
non-Gaussian CDM-dominated cosmological models. Our major conclusions include: 
(1) the epoch of reheating depends most strongly on the power spectrum of density 
fluctuations at small scales; (2) because of the effects of gas clumping, full reionization 
occurs at z ~ 10 in all models; (3) the CMBR polarization and the stars and quasars to 
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baryons ratio are strong potential discrimants between different assumed power spectra; 
(4) the formation of galactic spheroids may be regulated by the evolution of reheating 
through feedback, so that the Jeams mass tracks the non-linear mass scale; and (5) the 
evolution of the bias of luminous objects can potentially discriminate strongly between 
Gaussian and non-Gaussian PDFs. 

Subject headings: cosmology: theory — galaxies: formation — intergalactic medium - 
methods: analytical 

1. Introduction 

In this paper we model semi-analytically the thermal and ionization history of the universe 
and its effects on the formation of luminous objects. Many parts of this complex early evolution of 
the universe have been studied previously, beginning with the papers of Binney (1977), Silk (1977), 
and Rees and Ostriker (1977), the last of whom used models of the thermal history of the universe 
to investigate the formation of galaxies. Tegmark et al. (1997) used a simplified model of dark 
matter halo profiles to derive the minimum mass of objects which are able to cool in a Hubble 
time. The evolution of the intergalactic medium (IGM) has been explored semi-analytically by 
numerous authors, including Shapiro, Giroux, and Babul (1994), and recently Valageas and Silk 
(1999). Also recently, Ciardi et al. (1999) use a semi-analytic model for galaxy and star formation 
in conjunction with high resolution N-body simulations. The detailed evolution of the IGM was 
addressed in a full hydro dynamical simulation by Ostriker and Gnedin (1996) and Gnedin and 
Ostriker (1997) (hereafter G097). These studies explicitly allowed for clumping of the IGM, but 
was only able to simulate a small volume (a cube of side 2/i _1 Mpc). Ostriker and Gnedin (1996) 
found that a mass resolution of about 1O 4 M was required to resolve early epochs of reheating and 
reionization. 

The more general problem of galaxy (as opposed to dark matter halo) formation has also been 
addressed in full hydrodynamical simulations by many investigators, including Katz, Weinberg, and 
Hernquist (1996), G097, and recently by Pearce et al. (1999). Pearce et al. (1999) simulate a large 
volume, but with much courser mass resolution than G097, and thus was not meant to address 
the details of reheating and reionization. Kauffmann, Nusser, and Steinmetz (1997) developed a 
"hybrid" scheme, using a semi-analytic model of gas dynamics and star formation in conjunction 
with N-body simulations of the formation of dark matter halos. In most cases, including all the 
hydrodynamical studies cited here, only Gaussian CDM have been considered. 

Here, we attempt to incorporate much of the essential physics modeled by G097 into a semi- 
analytic model for the reheating and reionization of the universe. The work is not explicitly three 
dimensional. Rather, we treat the multi-phase medium in a statistical fashion and treat collapse 
in the approximation of the Press-Schechter method (Press and Schechter 1974). We consider the 
universe as an evolving two-phase medium — one phase (I) being neutral and the other phase 
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ionized (If). Of course we still must make many assumptions as to which physical processes are 
important and which may be neglected. The main assumption we make is that UV photons from an 
early generation of stars and quasars are the main source of energy for the reheating and reionization 
of the universe 

Our model also includes the following physical processes: (1) gravitational collapse and virial- 
ization, (2) non-equilibrium atomic and molecular physics of hydrogen, (3) heating and ionization 
by UV photons, (4) cooling, (5) simulated star/quasar formation, (6) clumping, and (7) feedback 
between star formation and the thermal and ionization evolution of the universe. All of these pro- 
cesses are calculated within the context of the standard evolution of the background cosmological 
model. In determining the evolution of structure formation, we explore several models for the prob- 
ability density function (PDF) and power spectrum of the initial cosmological density field. Once 
the background cosmology and normalization of the power spectrum are fixed, our semi-analytic 
model has only two free parameters related to the efficiency of star formation and UV photon gen- 
eration. These two parameters we fix by calibrating our results with G097 and by requiring each 
cosmological scenario to reproduce appoximately the observed ionizing photon intensity z = 4. The 
model we develop neglects the following processes: (1) spatial variations and correlations beyond 
the two phases, (2) non-virialization shock heating from gravitational collapse and supernova ex- 
plosions, (3) local optical depths to UV photons and the frequency evolution of the photon energy 
spectrum, and (4) helium, dust, and heavy elements in the IGM. Because these effects are likely 
to be important at late times, we only consider our models a reasonable approximation for z > 3. 

A word on motivation may be useful here. In almost all respects the full numerical hydro- 
dynamical simulations provide a better physical treatment than we are implementing here; for 
example, the neglected processes (l)-(4) not treated by the current work are included in G097. 
The two related ways that this treatment is to be preferred are that it is not limited by numerical 
resolution and that it allows, to lowest order, for a multi-phase (ionized/neutral) IGM. Given the 
overall relative strengths of the full numerical treatment, one can ask why bother to do anything 
less? The answer is threefold: 

1. First, residual doubts about numerical resolution always remain whatever tests are made, so 
it is important to see how robust the conclusions of the numerical work are when resolution 
limits are removed. 

2. The full numerical approach is so expensive that only a very small number of models can 
be explored. The result is that one is again concerned with the questions of robustness and 
generality of the results to variations of the adopted model. How sensitive are conclusions to 
the assumed nature of the perturbations (Gaussian/non-Gaussian), amplitude of the power 
spectrum, etc. Semianalytic methods allow one to explore the parameter space. 

3. The full numerical results, paradoxically, in so far as they are better and are more and 
more comprehensive in physical treatment, tell one less and less! Understanding depends on 
knowing which physical processes and initial conditions are essential to the results, and which 
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merely modify the details of the overall picture. A semi-analytic treatment allows one to add 
and subtract from the modeling to determine what matters and what does not. 

The organization of this paper is as follows. In § 2, we present the elements of our adopted 
physical model for reheating and reionization. In § 3, we apply it to several Gaussian and non- 
Gaussian cosmological models. We then place some constraints on these models using several 
current and potential observations, such as Gunn-Peterson absorption and the CMBR, which are 
directly affected by the ionization history of the universe. In § 4, we discuss the relationship between 
the reheating and ionization history of different models to the evolution of the luminous matter 
in the universe. We summarize our results in § 5. In the appendices, we delineate our treatment 
of the standard physics of the universe, as well as the specific calculational implementation (e.g., 
finite difference equations) of our semi-analytic model. 



In this section, we describe the semi-analytic model we use of the physics of reionization. 
We concentrate on the treatment of a two-phase medium with sources and sinks of radiation. In 
particular we define our basic differential equations, and describe how we calculate the formation 
rate of bound objects, the UV production rate, and the clumping of cosmic gas. 



Since we are considering a two-phase medium, we will need to understand the time evolution 
of the filling factor of ionized regions ///. The two other independent quantities required are the 
ionization fraction in the ionized region xu and the mean (volume-averaged) UV ionizing (> 13.6 eV 
= eo) intensity J. The corresponding intensity in the ionized region is Jn = J/ fn- in our notation, 
Tin is the total mean hydrogen density. 

Considering only ionization of hydrogen, local particle conservation in the ionized region gives 



2. A Model for the Physics of Reionization 



2.1. Definitions and Conservation Equations 



^-(n H xii) = -3-n H xn + n H (l - xii) ^ Ju<Jv 
at a en 




(1) 



(2) 



where the first term on the right hand side is from the expansion of the universe, and the second 
term from photoionization, and the third term from collisional ionizations (rate coefficient k c i) and 
radiative recombinations (rate coefficient ajj). The effective cross section to photoionization a p is 
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given below. The clumping factor C = (pf) / (pb) 2 takes into account the non-uniformity of the gas. 
We describe how to estimate C in § 2.4. 

For an energy spectrum of the form E v oc z^~ 7 , the average UV energy density (E = 4ttJ/c) of 
the universe evolves as 

^ = 5-[4 + ( 7 -l)]-£- fnn H (l-x u )47rJ n a e -eon H x u ^- (3) 
at a at 

CL df 
= S - [4 + (7 - 1)]-E ~ «h(1 - x n )Eca e - eQn H x u — , (4) 

where we have used the relationship between E, Ejj, and Jn- We use a value of 7 = 1.5 throughout. 
The integrated source function is given by S, the Hubble expansion is accounted for in the second 
term, energy lost to maintaining the ionization of the ionized region is given by the third term, and 
the last term is due to energy lost to ionizing new (i.e., previously neutral) material. The effective 
UV energy loss cross section is a e . 

The definitions of a p and a e are as follows. If the ionizing specific intensity is J u , then the 
ionizing intensity is 

J = ( * Judv, (5) 



J l/ 

where is a high frequency cutoff (we use 7.4 keV, from Abel 1998). For a power law 

J, = «W(v) 7 ' (6) 

we have the relation 



1 - 



V * 



7-1 

The photoionization rate per hydrogen atom is given (in s _1 ) by 



(7) 



T p = r ^3 v o{v)^- (8) 

_ Mo j- my * 



where h p is Planck's constant. Equating this with 4ir Ja p /eo gives 

Similarly, assuming that all the excess energy from each photoionization goes into heat and 
produces no UV photons, the rate of UV energy loss per hydrogen atom (erg s" 1 ) is 



r e = 



= / 4nJ v a(v)dv (11) 

= 4,j rr-°y*^. (12) 
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Equating this with Air Ja e gives 




The difference between these cross sections gives the cross section for photoheating, or equiv- 
alently a photoheating rate (erg s _1 ) of 

T ph = 4irJ(a e -a p ). (14) 

In order to fully specify the model, we need to be able to determine S, as well as have a 
relationship between E and ///. These turn out to be intimately related. 



2.2. Formation Rate of Bound Objects 

We assume that stars/quasars only form in virialized collapsed objects. Hence the first input 
into the star formation rate is the formation rate of bound objects. In the next subsection we relate 
these bound objects to UV energy output in order to derive expressions for the source function S 
and the relationship between the mean radiation energy density E and the filling factor fjj. 

Generalized Press-Schechter for a PDF P(v), and normalization factor f^ 1 = P{v)dv 
gives the comoving number density of collapsed objects with initial comoving radius between R 
and R + dR to be 

» T ,^ n „r m 3 dlnv(R,t) ,^.dR , . 

N R iR = f .P H R, t)] —-^ ]t > v{R , t) - (15) 

- (R ' () s ^kJW) (16) 

dlnv(R,t) d\n<T-\R) _ 

ainii = rflnB (17) 

where a(R) is the rms fluctuation on scale R linearly extrapolated to the present, S c ps 1.69 is 
the standard linear overdensity of collapse, and D(t) is the linear growth factor normalized to 
D(to) = 1. The original Press-Schechter (Press and Schechter 1974) method was derived in for a 
Gaussian PDF, and has been checked in that case with N-body simulations (Efstathiou et al. 1988). 
These expressions in the non-Gaussian case were derived in Chiu, Ostriker, and Strauss (1997), and 
recently tested through N-body simulations for several non-Gaussian models (Robinson and Baker 
1999). 

For each model, we filter the power spectrum using a top-hat window to obtain the rms 
fluctuation a. The power spectra we use are from Bardeen et al. (1986) and Sugiyama (1995) in 
the adiabatic cases; Pen, Spergel, and Turok (1994) and Pen and Spergel (1995) in the textures 
case; and Peebles 1999ab in the ICDM case. 
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In general, we cannot use only the Press-Schechter number density to determine the formation 
rate of bound objects. This is because the net time-rate of change of the number density of objects 
at a particular scale consists of both a formation Nrjotui rate and a destruction rate NR^est- Thus 
the actual formation rate may be larger than the rate inferred from Press-Schecter alone. Of course, 
we must constrain these rates to be consistent with the Press-Schechter formalism: 

Nr = Nrm™ ~ ^,dest, (18) 

where Nr is given by differentiating equations (15)-(17) with respect to time. We seek a semi- 
analytic methodology applicable to an arbitrary PDF which gives these two rates. 

Sasaki (1994) models the formation and destruction rates of these objects by assuming that 
the destruction efficiency rate 4> = A^dest/-^ has no characteristic scale — "destruction" by 
incorporation into larger objects is hierarchical and scale-invariant. In this section, we generalize 
his derivation to non-Gaussian PDFs. 

The assumption of no characteristic scale of destruction efficiency, in addition to the require- 
ment that the net rate of formation minus destruction be consistent with Press-Schechter, leads to 
the conclusion that (f> must be a function only of time. Sasaki (1994) derives the following relation 
for the destruction probability per unit time: 

m = (19) 

which we find to be independent of the PDF. It follows that the probability that an object created 
at t± exists at ti without having been incorporated into a larger object is given by 



Psurv(*2|*l) = exp 



rt2 

/ 4>(t)dt 
Jtl 



m 



Combining 4> and the Press-Schechter consistency requirement equation (18) leads to a forma- 
tion rate of the form: 

D f-P'\ 

N RMm = N R ■ u(R, t)- [— , (21) 



D V P 

where P' = dPjdv. The quantity —P'/P for a Gaussian and an exponential P{v) takes on the 
forms: 

~ P v, if P(u) oc e-^ 2 ; (22) 



P 

= u, if P{u) oc e- uv . (23) 

Textures and non-Gaussian ICDM both have PDFs which are well approximated by exponentials 
(Gooding et al. 1992; Park, Spergel, & Turok 1991; Peebles 1998ab). For textures u = 1.45 (using 
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Park, Spergel, & Turok 1991 as described in Chiu, Ostriker, & Strauss 1997), and for non-Gaussian 
ICDM uj = 0.67 Peebles (1999b). Recalling that 

§-'<»>;• < 24 > 

where the velocity factor f(fl) ~ SI ' 6 , the formation rate in closed form is 

NRMmdR = f- [-P»] ^/(^^(^^(^t)^. (25) 

This is the formation rate per comoving volume we use in our model. For a power law a(R) oc R~ a , 
the logarithmic slope d\nv/d\nR = a is a constant. 

This formalism for merging and formation rates is consistent with the work of Blain and Lon- 
gair (1993), who used a numerical extension of Press-Schechter with simple forms for the merging 
probability. All other work of this type is based on Lacey and Cole (1993) and Lacey and Cole 
(1994), who construct conditional probabilities by extending the work of Bond et al. (1991). Un- 
fortunately, this method is heavily dependent upon the PDF begin Gaussian, and is sufficiently 
complicated so that extending it to non-Gaussian PDF is not straight forward. Furthermore, since 
objects are "constantly growing" in the Lacey and Cole (1993) formalism, rather than in "punctu- 
ated equilibrium" as implied by equation (20), it is difficult to directly compare the two methods. 
However, one sign that they might be roughly consistent is that Viana and Liddle (1996) obtain 
similar results for cluster abundances whether they use the Lacey and Cole (1993) or the Sasaki 
(1994) methodology. 

We further note that during reheating and reionization, structure formation is largely hier- 
archical in most of the models we consider. This is because only a small fraction of the universe 
needs to collapse in order to release enough energy in order to reheat and reionize the universe. 
Collapsed objects are thus still "rare" (v > 1), and their net formation rates are thus dominated 
by their creation rates. The use of the Sasaki results thus provide only a small correction. In the 
limit where the "destruction" rates are negligible, all semi-analytic merging formalisms which are 
based on Press-Schechter must give the same answer in order to be consistent with equation (18). 

2.3. Determining the UV Production Rate 

Our basic model for UV production is that individual halos which can undergo star forma- 
tion will emit UV radiation radially, and thus be surrounded by a sphere of ionized material, a 
"cosmological Stromgren sphere." Consider halos of baryonic mass in the range M& to + dMf,, 
corresponding to comoving length scale between R and R + dR, formed in the time interval t to 
t + dt, and "observed" at time t. The comoving number density of such objects is given by 



N(R, t, t) dR dt = N iorm (R, t) ■ p SU rv(*|*) dR dt, 



(26) 
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where the notation is the same as in the previous section. The second factor accounts for the 
destruction objects of scale R through merging into larger objects. Although we are using Sasaki's 
method to determine -/Vf orm and p surv as described above, the form of equation (26) is completely 
general. Note that Mb and R are related by Mb = AttR 3 pbfi/3. Following the prescriptions of Cen 
(1992) and Gnedin (1996), we assume that star formation is spread over a dynamical time scale, 
so that the UV luminosity of such an object is given by 



i(M b ,t,t) = (i - /,)M 6 c 2 £cooie ; ehmeuv (i-J. ) , X p 

Myn V Myn 



t-t 



yn 



(27) 



where /* is the fraction of the baryonic mass which has already been processed into stars/quasars 
and tdyn is the dynamical time associated with the object at its time of formation (see equation 
CI). We are also assuming here that the stars/quasars are uniformly distributed throughout the 
collapsed halo. 

The efficiencies e^ are as follows: e coo i is the mass fraction of the halo which is cooling, e* is 
the mass fraction of cooling gas which physically ends up in collapsed objects, eh m is the mass 
fraction of collapsed objects which are quasars or high-mass stars (i.e., with substantial UV flux), 
and euv is the mass-to- UV efficiency of the high-mass stars/quasars. The cooling fraction e coo i we 
derive below in Appendix C, and is a function of the halo virial mass and the redshift of collapse. 
G097 uses euv = 6 x 10~ 5 and Gnedin (1996) gives eh m = 0.16, derived using a standard initial 
mass function. We cannot use exjv directly because we are not resolving the mass elements within 
individual cooling halos. The "resolution" factor e* parameterizes this uncertainty. We use a net 
efficiency 

e e ff = e*eh m euV; (28) 

normalized by requiring, for each model, that the ionizing intensity J21 = 1 at z = 4. There is, 
effectively, only one free parameter here, e e ff , and that is fixed by normalizing to observations. 

A separate estimate of e*, however, is necessary in order to derive /*, the fraction of baryons 
which is sequestered in stars/quasars and their remnants. It is important to keep track of /* to avoid 
"double counting" star formation — baryons which are in such bound objects are not available for 
additional star formation. Note that results are insensitive to the value of /* while /* < 10%, since 
they are only dependent on 1 — We can use the G097 simulation to derive our "resolution" 
factor e* by assuming that it is constant for all cosmological models. If e* is a constant, then it 
should correspond to the value which gives the same J21 as the G097 simulation, using their values 
for ehm and euv and the same cosmological model. This correspondence leads to a value e* ps 0.5. 

Given e*, the simplest way to calculate the fraction /* is actually to integrate over the UV 
production rate S, given in equation (31), dividing by the appropriate efficiencies: 

rt 



m = — (29) 

r * S{t)dt 



I 

Jo 



e eS nHm H c 2 / 'X' 



(30) 
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where X is the hydrogen mass fraction, and we have expressed the result in terms of the J21- 
normalized effective efficiency e e g- and the resolution factor e*. As a check on our value of e*, we 
compare /* at z ~ 4 between our model and the G097 simulation, given the same cosmological 
parameters. We find the values to be consistent, with /* ~ 2% at this redshift (see also Figure 2). 

The global mean UV production rate per unit (physical) volume is given by 

/'OO ft 

S(t) = a(ty 3 / dR dtiP(R, t) ■ N(R, t, t) ■ £(M b , t, t), (31) 
Jo Jo 

where the function ip(R,t) is a selection function for objects which are able to collapse and form 
stars/quasars. One simple form for tp is 

Aim P {R,t) = 0[R-R,j{t)}, (32) 

where the 9 function is zero if the argument is < and unity if the argument is > 0. These represent 
the criteria that a mass must be greater than the Jeans mass to collapse. In order to account for 
the fact that the Jeans mass will be different in the two phases, we can define ip to be 

V>2 P h(#, t) = { (1 - ~f u ) ■ 0[R — Rj, c (t)] +f n -e[R- Rj, h (t)] } , (33) 

where fn is the filling factor of the ionized region at time i, Mj jC is the Jeans mass in the neutral 
(cold) region, and Mj^ is the Jeans mass in the ionized (hot) region. Note that we are assuming 
that an object, once formed, can only be destroyed by merging and forming a larger object. We 
are assuming, thus, that virialized gas halos are "self-shielded" from ionizing radiation due to their 
high density. That is, they can not be heated to the point of dynamical instability nor can they be 
"evaporated" away. 

In order to derive the needed relation between fn(t) and E(t), we make two further assump- 
tions. One is that the ionized volume vn(Mb,t,t) around an isolated object is proportional to its 
UV luminosity £(Mb, t, t) with a coefficient of proportionality n(t) that is a function only of the 
time t when it is "observed." This should be a reasonable approximation, since the recombination 
rates are determined by the density and temperature at t, not at the formation time t. As long 
as objects live longer than the recombination time, our assumption is consistent with the work of 
Shapiro and Giroux (1987) on cosmological H II regions. We thus have the relation 

£(M b ,i,t) = K-\t)v n (M b ,t,t) = 3^y4- (34) 

The second assumption, which is needed to determine K.(t) from ///(t) and S(t), is that the 
filling factor can be related to the average number of "overlapping" ionized regions by modeling 
it as an uncorrelated Poisson process. The average number X(t) of overlapping regions seen by a 
random volume element is 

f'OO ft 

X(t)= dR dtip(R,t) ■ N(R,i,t) ■ v n (M b ,t,t). (35) 
Jo Jo 
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We no longer have the factor a 3 because vn is a comoving volume. For Poisson probabilities, the 
filling factor fn(t) is related to X(t) by 



f II (t) = l-e^ t l 
Combining equation (35) with equation (34) and using (31) gives 

\(t) = n(t) a 3 {t)S(t). 
Thus we have determined n(t) in terms of fnit) and S(t). 



(36) 



(37) 



With K,(t) defined, we can derive a relationship between fn(t) and E(t) which will close our 
system of equations. Consider the contribution of each source to the mean intensity J(t). Neglecting 
the optical depth within an isolated ionized region, and assuming no flux escapes beyond a radius 
rn(M b ,t,t), the average flux seen by a random volume element due to an object of luminosity 
£(M b ,t,t) is 



d[4vrJ(t)] 



dRdtN(R,t,t) [" d 3 r «nM) 
Jo 



4vr[a(t)r] 2 

dRdtN(R, t, t) ■ £(M b , t, t) ■ r n (M b , t, t) ■ a~ 2 (t) 



dR dtN(R, t, t) ■ £* /3 (M b , t, t) 



Air 



nl/3 



a- 2 (t), 



where r is a comoving length, a r is a physical length, and we have used 



r n (M b ,t,t) = 



3n(t) 

47T 



£(M b ,i,t) 



nl/3 



(38) 
(39) 
(40) 

(41) 



which we assume to be much smaller than the Hubble radius. Equation (40) also assumes that 
the central luminous object emits radiation radially, which is a good approximation as long as the 
radius of the sphere is much larger than the virial radius. Using equation (37) in equation (40) 
gives 



d(AitJ) = dR diN{R, t, t) ■ £ 4/3 (M b , t, t) 



3A(t) "1 1/3 



a~ 3 (t). 



.4vrS(i)_ 

Therefore the mean energy density E(t) = 4irJ(t)/c at time t is given by 

1/3 



(42) 



E(t) = - 



3A(t) 
c [4irS{t) 



1 / dR dti/j(R, t) ■ N(R, t, t) ■ £ i/3 (M b , t, t). (43) 
Jo Jo 



This is the desired relationship between fn(t) = (1 — e A( ^) and E(t) which closes the system of 
equations. 

We have neglected spatial correlations among these sources. This approximation should be 
less important for the texture model, which is well approximated by uncorrelated seeds, than for 
the Gaussian models. As long as the filling factor is small, though, these correlations are likely not 
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to be important, since most of the universe is neutral. We will find below that in the Gaussian 
cases we consider, reionization occurs quite rapidly, and this rapid percolation likely overwhelms 
any spatial correlations, which tend to slow down the percolation. 

Furthermore, we have neglected the optical depths within the collapsed object as well as in 
the gas in the Stromgren sphere surrounding each object. We assume both that the covering factor 
of other halos within the sphere is small, and that there is a sharp transition between the ionized 
and neutral regions. This approximation is somewhat alleviated by assuming a clumping factor for 
the gas (see § 2.4). With clumping, our model is self-consistent in the volume average, accounting 
for all photoionizations and recombinations in the ionized region. 

In summary, we assume that star formation in collapsed halos is the source of UV photons. 
The number density of these sources is determined using the Press-Schechter formalism with the 
requirement that the baryonic halo mass is greater than the Jeans mass; the luminosity of each 
source is determined by the fraction of cooling gas and several efficiency factors. These sources 
are assumed to be scattered randomly (i.e., a Poisson point process), and each surrounded by a 
Stromgren sphere. The volume emmissivity of these objects is assumed to be uniform, and they 
are assumed to radiate radially into the surrounding intergalactic medium. Assuming that at fixed 
time, the volume of each Stromgren sphere is proportional to the luminosity of the central source, 
the radius of each sphere is determined by requiring consistency among the filling factor, the mean 
intensity, the number density of sources, and the total UV luminosity. 

2.4. Nonuniformity of Cosmic Gas: The Clumping Factor and Tracking the Diffuse 



The importance of clumping in determining the ionization history of the universe was stressed 
by G097. Because the clumping factor has been neglected in most previous semi-analytic work 
on reionization, there has been no attempt to estimate it analytically. In this subsection we use 
Press-Schechter to estimate the clumping factor. A similar procedure was described in Valageas 
and Silk (1999). 

Consider dividing the baryons in the universe into clumped and diffuse components. The 
clumped component consists of baryons which have collapsed into virialized halos. Using equation 
(26), the fraction of the baryonic mass of the universe which collapsed df m in time interval t to 
t + dt, and survives without being incorporated into a larger object to time t, is given by 



where ip(R,t) is the selection function which filters out objects below the Jeans scale. The fraction 
of baryonic mass which is in collapsed objects at time t is given by integrating over t: 



Component 




(44) 




(45) 
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Note that we do not take the collapsed fraction f m (t) to be equal to be the fraction of mass above 
the Jeans mass at time t. This is because during reheating, the Jeans mass increases, but objects 
below the Jeans mass which formed at earlier times (when the Jeans mass was lower) can still 
survive as long as they are not merged into larger objects. If all of these objects are in virialized 
halos, then their overdensity relative to the mean density will by 

Pvir _ r _ ^vir 

— = 6 ™ = wr m 

where A vir is the overdensity relative to the critical density. Thus the fraction of the volume in the 
universe that they occupy will be 

/„(*) = j^- (47) 

O v ir 

Self-consistency requires that the diffuse component be underdense so that pais < P- Taking volume 
average gives 

P = fvpvii + (1 - fv)Pdis, (48) 



which implies that 

Pdiff 1 - /r? 



(49) 



P 2 [ fmKir + ( \ f 7 )■ (51) 



P 1 - fv 

The volume-averaged square density will be 

(P 2 ) = fvP 2 ^ + (l-fv)p 2 diS (50) 

1 fm 
l-fv 

The average clumping factor C is thus 

C ee (p 2 )/p 2 (52) 
= fmS^+ ^'J^ 2 . (53) 

This is the clumping factor we use in our model. Our model, however, assumes that UV sources 
are uniformly distributed throughout collapsed halos and that all halos carry the same overdensity. 
These two assumptions will tend to reduce clumping, so equation (53) probably gives a lower bound 
on the importance of clumping. As shown in Figure 2, the results of this semi-analytic calculation 
are in fair agreement with the results of the numerical simulation in G097. 

We keep track of the ionization and temperature history of the purely diffuse phase of the 
ionized region. We do this in order to properly calculate the Gunn-Peterson optical depth, which 
probes the diffuse intergalactic medium (IGM). This phase is by definition unclumped, and is 
underdense relative to the mean, equation (49). Note, however, that this diffuse phase plays no 
feedback role — it simply responds to the other physics which is going on. 
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3. Results 

In this section we present the results of our calculations for a range of cosmological models. 
We begin with a "calibration" model which uses the same cosmological parameters as were used by 
G097 (fio = 0.35, Q A = 0.65, h = 0.70, and Sl b = 0.03). We then consider adiabatic Gaussian CDM 
(G+ACDM), textures (T+CDM), and non-Gaussian isocurvature CDM (ICDM) as proposed by- 
Peebles (1999ab). To help assess the sensitivity of the results to the PDF and the power spectrum, 
we also consider an "artificial" case of a texture PDF with an adiabatic CDM power spectrum 
(T+ACDM). We use a Hubble parameter h = 0.65, a baryon density parameter Sl^h 2 = 0.0125 
throughout, and consider flat (A-dominated) and open models with matter density parameters of 
S7o = 1, 0.45, 0.35, and 0.25. All models are normalized to cluster abundances (method of Chiu, 
Ostriker, & Strauss 1997 for all models except ICDM, which uses Peebles 1999ab). The cosmological 
parameters associated with the models we consider are listed in the first three columns of Table 1. 
As noted above, the UV efficiency e e ff in each run is adjusted to reproduce J21 = 1 at z = 4. There 
are no other adjustable parameters in the model. 

3.1. Comparison with G097 Simulation 

In Figure 1 we show the filling factor ///, the ionizing intensity J21, the neutral fraction (1 — x), 
and the gas temperature T in the case were we mimic the G097 calculation. While the ionizing 
intensity was adjusted to be consistent with the G097 value at z = 4, the evolution of J21 is similar 
in the two calculations. As in the G097 calculation, reionization to 1 — x = 10 -3 occurs suddenly 
at z ~ 10, while reheating occurs more gradually. The main difference between this work and G097 
is that reionization and reheating begin somewhat earlier in our model. This may in part be due 
to the limited resolution and single-phase treatment of G097. 

In Figure 2, we show the stellar fraction and clumping factor for this reionization calculation. 
Again, the values are similar to those in G097. The clumping factor shows the greatest discrepancy, 
and is without doubt inaccurately calculated in our model. Without full-scale simulation, however, 
it is difficult to calculate the clumping factor to a much greater accuracy than what we have 
presented in § 2.4. The fact that our result is within a factor of a few of the G097 calculation is 
encouraging. The earlier onset of star formation shown in the top panel by our semi-analytic model 
may reflect the limited resolution of G097, which while adequate for average objects, is insufficient 
for early, dense small ones. 

In Figure 3, we show the nonlinear and average Jeans mass scales derived from our model. Our 
results are qualitatively consistent with OG96: we find that the Jeans mass and the nonlinear mass 
(here defined by 5M/M « 0.4 ^ track each other to remarkable accuracy during the reheating phase. 
Because reheating begins earlier in this model than in the G097 simulation, the precise degree of 
nonlinearity which tracks the Jeans mass is different, and the epochs during which the Jeans and 
nonlinear mass scales track each other occur earlier. A possible physical basis for this tracking is 
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clear. When the Jeans mass is lower than the nonlinear mass, star formation accelerates so that 
the average temperature increases, the Jeans mass increases to above the nonlinear mass scale, and 
further star formation is quenched. 

We do not expect an exact correspondence between our method and the G097 work. For 
example, this semi-analytic calculation is not limited by mass resolution, so mass condensations 
down to the Jeans mass limit can contribute to the ionizing UV flux. This may partially explain 
why reheating starts at z ~ 30 in our model, rather than at z ~ 20 in G097. Another contributing 
factor to this difference in the reheating epoch is that the largest G097 simulation lacked significant 
large-scale power relative to the true power spectrum. Furthermore, we model here a two-phase 
medium where the radiation field is not spatially uniform, while G097 use an average UV radiation 
field. On the other hand, this calculation neglects all spatial correlations, which are automatically 
taken into account by the three-dimensional numerical calculation of G097, and uses a power-law 
radiation spectrum. While it is unclear what the net effect of all these differences should be, it is 
encouraging that two completely different methods of calculating the reheating and reionization of 
the universe give such similar results. 



3.2. Reheating and Reionization in Various Cosmologies 

The results of our calculations for the models we consider are presented in Table 1 and Fig- 
ures 4-9. Qualitatively, the evolution of reheating and reionization can be divided into two broad 
categories characterized by late or early reheating, corresponding to late or early evolution of 
the filling factor fjj. In models with late reheating, which include the Gaussian models and the 
T+ACDM model, the filling factor rises from f n < 10~ 4 at z = 30 ~ 50 to f n ~ 1 at z = 10 ~ 20. 
For models with early reheating, which include the pure texture T+CDM models and the ICDM 
model, the filling factor is already fn > 0.1 at z ~ 80, but rises slowly to fn = 1 at z = 10 ~ 20. 

The general picture of the evolution in these models is as follows. We note that by "reioniza- 
tion" we mean when fjj is approaching unity, but that by "full reionization" we mean the epoch at 
which the individual ionized regions start overlapping substantially, so that A = — In (1 — fn) 3> 1. 

For an adiabatic CDM-type power spectrum, the fluctuation spectrum ctr rises logarithmically, 
or nearly logarithmically, as R — > 0. Hence, at the Jeans scale of ~ 0.017i _1 Mpc, there is a 
logarithmic cutoff in power for z > (JRj ~ 10 — 20. Furthermore, cooling becomes less efficient at 
the Jeans scale at later times. These two factors lead to a suppression of reheating and reionization 
until late times. Full reionization occurs after the filling factor approaches unity, at redshifts z < 10. 

These results for Gaussian CDM models are consistent with the work of G097 (who considered 
a flat Qa = 0.35 CDM model), as well as Ciardi et al. (1999) (who considered an Qo = 1 CDM 
model), but vary somewhat from the results of Valageas and Silk (1999) (who considered Qq = 1 and 
open Slo = 0.3 CDM models). The latter find a slightly lower value for the redshift of reionization 
(z ~ 6 — 7). This difference is unlikely due to gas clumping, as proposed by Ciardi et al. (1999), since 
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both this work and G097 incorporate clumping. In fact, our work uses a prescription for clumping 
very similar to Valageas and Silk (1999). We do, however, use a different prescription for calculating 
the radius of the cosmological Stromgren spheres, but the rapidity with which reionization occurs 
may mitigate the effect of this differences. While the exact value of the epoch of reionization 
remains uncertain, various methods, including ours, give the range z = 8 ± 2 for Gaussian CDM 
models. 

For T+CDM and ICDM power spectra, the fluctuation spectrum rises as a power law for small 
R, the Jeans scale becoming nonlinear with SM/M ~ 1 at z ~ 50. At this redshift, cooling at the 
Jeans scale is still efficient, and UV photons are created in abundance. At z ~ 100, the Jeans scale 
is already mildly nonlinear, with 5M/M ~ 0.5, causing the universe to already be substantially 
reionized (/// > 10%) at z ~ 80. But because of the early formation of non-linear structures, the 
clumping factor is higher, and hence the recombination rate is greater. Full reionization is thus 
delayed until z = 10 ~ 20. Detailed calculations for reheating and reionization have not been done 
for these models, so comparisons with other work can not be made. 



3.3. Current and Future Constraints on the CMBR: The Compton y and the 

Optical Depth to Recombination 

The y-parameter measures the CMBR spectral distortion away from a blackbody spectrum 
caused by Compton-Thomson scattering. Free electrons scatter and so Doppler shift CMBR pho- 
tons, causing them to undertake a random walk in frequency. In thermal equilibrium, of course, the 
spectrum does not shift, since by detailed balance, all the processes must cancel. However, after 
thermal decoupling at z ~ 100, the kinetic temperature of the gas evolves away from the CMBR 
temperature, first cooling by adiabatic expansion, and then heating by virializing in collaping 
structures or by photoheating from star formation. The parameter y is defined by: 

• fc-B (T e — T) 

V = n °Tn e c, (54) 

m e cr 

where T e is the electron temperature, T is the CMBR temperature, and the other symbols have 
their usual meanings. When structure formation begins in earnest and the gas temperature rises, 
low-energy CMBR photons will be scattered to higher energies. 

The y we calculate in our models using equation (54) are shown in Table 1. Again, the T+CDM 
and ICDM models behave similarly, with y = 2 ~ 11 x 10~ 7 , and the Gaussian and T+ACDM 
models have y = 1 ~ 3 x 10~ 7 . Thus, they differ by a factor of ~ 3. The COBE satellite (Fixsen 
et al. 1997) put a bound on the spectral distortion 

y < 1.5 x 10" 5 . (55) 



At long- wavelengths, (in the Rayleigh-Jeans limit) the effective temperature T e g = T+5T is lowered 
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by a factor proportional to y: 

Y = -2y. (56) 

Unfortunately, the prospect for future experimental bounds on y which are significantly better than 
equation (55) are not good (Wilkinson, D., 1998, private communication). 

As a check, an upper bound on y can be obtained by assuming that all the energy released by 
star formation before redshift z is transfered to the CMBR. If a fraction de of the baryon rest mass 
which has been transfered, then the contribution the y-parameter, dy, is governed by 

dep b = idy p 1 , (57) 

or 

dy = de^- (58) 

4/9 7 

* ^ : 10-4(1 + z) (59) 

where p 1 is the photon energy (mass) density 

Compton cooling is only efficient (i.e., faster than the expansion rate) when z > 5. For the 
Gaussian models, the mass fraction in stars/quasars at this epoch is /* ~ 0.03, and using the UV 
efficiency euv ~ 6 x 10 -5 gives a bound e < 2 x 10 -6 , which implies that y < 4 x 10 -5 . This 
assumes that all the energy is transferred at z ~ 5. Integrating equation (59) gives a value of 
y < 0.6 - 1.3 x 10" 5 , all of them below the COBE limit. The numbers for the T+CDM and ICDM 
models are about y < 0.3 — 0.7 x 10 -5 . The remaining factor of ~ xq 16-2 ' 6 comes from several 
factors. One is the fact that the ionizing photons leave in the electrons on average ~ 1 — 2 eV 
per ionization in heat which can be transferred to the CMBR, resulting in a ~ 10% efficiency 
factor. Furthermore, there is competition from cooling processes other than Compton cooling, such 
as H2 and atomic line cooling. Thus the CMBR only absorbs a small fraction of the total UV 
energy produced. It appears, then, that the y-parameter is not a good distinguishing characteristic 
between these models. 

A potentially stronger constraint on reionization comes from the CMBR anisotropy spectrum. 
The optical depth to recombination is given by 

n e caTdt (60) 

On small angular scales, reionization results in the damping of the temperature anisotropy spectrum 
by a factor e~ Tlcc . Because we stop our calculation at z = 3, we can only give a lower bound on 
the recombination optical depth. 

Zaldarriaga (1998) notes that while it might be difficult to extract the difference between a 
small degree of damping due to reionization and changes in other cosmological parameters, even 
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a modest optical depth due to reionization can enhance the polarization anisotropy spectrum on 
medium angular scales by an order of magnitude. Reionization causes a polarization anisotropy 
because the polarization source is dominated by the temperature quadrupole, which grows after 
recombination via free-streaming of the temperature monopole. For example, Zaldarriaga (1998) 
shows that an optical depth to recombination ~ 0.5 can produce a polarization anisotropy at 
I = 10 on the order of a micro-Kelvin (/xK), while the same model with T rec = gives a polarization 
anisotropy of of ~ 10 _L5 //K. 

In the models we consider, those with the adiabatic CDM power spectrum give a small optical 
depth of r rcc < 0.1, while the other models have T rec > 0.3. According to Zaldarriaga, a negative 
detection on the 1° scale with sensitivity ~ 1//K would be able to rule out these moderate optical 
depth models. The upcoming MAP and PLANCK satellites may be able to distinguish between 
the low and high reionziation optical depths. As long as polarization foregrounds do not limit 
the measurements, the MAP and PLANCK satellites should be able to detect reionization optical 
depths of ~ 0.02 and 0.004, respectively (Spergel, D., private communication). 



3.4. Constraints on the Intergalactic Medium: The Gunn-Peterson Optical Depth 

and the Baryon Fraction 

Gunn-Peterson optical depth as a function of scalefactor a is given by (Peebles 1993) 

3A 2 iA| 1 n H o (a)c 

T~GP = , 61) 

AH a V^V~ 3 + + (1 - - 0\)a" 2 

where A21 = 6.25 x 10 s s -1 is the Lya decay rate and A21 = 1216 A is the wavelength. The Gunn- 
Peterson optical depth at z = 3 is shown in Table 1. The results for all models are near the upper 
limits of 0.1 ~ 0.35 at z = 4 established by observations (Jenkins and Ostriker 1991; Giallongo et 
al. 1994). Because we have neglected the effect of large-scale shock heating (which would increase 
the temperature) as well as metals and other coolants (which would increase the star formation 
rate), it is plausible that our final neutral fractions are overestimated, leading to an overestimate 
of r G p. 

A stronger discriminant between models may be in their predicted star formation histories. 
Figure 10 shows the baryon fraction in stars/quasars predicted by the various models we have 
considered. Clearly, the early reheating models require that a significantly greater fraction of the 
baryons be tied up in stars/quasars or their remnants. The observational limit that we consider 
is from the analysis by Rauch et al. (1997). Their work indicates that Lya clouds account for the 
majority of the baryons in the universe. The strong limit on the Lya clouds of Q c \h 2 > 0.017 
combined with the most generous (low deuterium abundance) big bang nucleosynthesis (BBN) 
constraint value of f^/i 2 = 0.024, implies that a fraction less than 0.017/0.024 = 0.29 of the 
baryons in the universe can have been converted into stars/quasars by a redshift z = 3. The 
models with the adiabatic power spectrum are well within this constraint, with stellar fractions 
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/* < 0.06. The T+ACDM models gives lower values than the Gaussian models because of the 
lower value of cluster normalization required. The ICDM model is within the limit with /* = 0.2; 
however, if the BBN value of Qbh? falls below ~ 0.02, or if the limit is increased by more than 
~ 15%, then the ICDM model will also encounter difficulties. The pure texture T+CDM gives 
stellar fractions which are at or above the limit set by Rauch et al. (1997). 

While we have used Vt^h 2 = 0.0125 in these models, we expect the fraction of baryons converted 
to stars/quasars should be similar or higher for models with higher Q^- The stellar fraction depends 
on the fraction of baryons collapsing and on the fraction of collapsed mass which is cooling. In- 
creasing the baryon density will increase the cooling efficiency, thus decreasing the Jeans mass and 
increasing the cooling fraction. Both of these effects will lead to an increased fraction of baryons un- 
dergoing star formation and thus an increased stellar fraction. Of course, the mass-to-UV efficiency 
will need to decrease so as to produce the same J21 at z = 4. 

4. Discussion: Physical Bias and the Formation of Galactic Spheroids 

In this section, we explore the relationship between luminous matter and the thermal and 
ionization history of the universe. We begin by addressing the issue of physical bias at the Jeans 
mass limit using a generalized form of statistical bias based on the peak-background split. We then 
consider the cooling efficiency of objects at or above the Jeans mass, and speculate as to the origin 
of galactic spheroids. 



The peak-background split was first suggested by Kaiser (1984) and recently revisited in the 
case of non-Gaussian models by Robinson, Gawiser, and Silk (1998). In this model, the density 
contrast field 5{x) is split into a short-wavelength "signal" S s and a long-wavelength "background" 
5b- The split is made in such a manner that S s and Sj, are largely uncorrelated, (5<A) « 0. The 
wavelengths of 8 S corresponds to the scales of the objects whose clustering properties are of interest, 
and the wavelengths of 5{, corresponds to the scales at which the clustering properties (and bias) 
are to be measured. In this prescription, the local threshold for collapse S c is modified locally to 
<5c — <5f>( x )) and the (observed) Eulerian density is modified from the Lagrangian one by a factor 
1 + <5&(x). Let the average comoving number density of objects of scale R to R + dR be given by 
Nr(5 c )cIR, where we have explicitly shown the dependence on the critical threshold of collapse S c . 
Then the local comoving number density per unit scale dR is given by 



4.1. 



Generalized Statistical Bias 



[1 + 4(x)] • N R [6 c - <5 6 (x)] 



(62) 



= [1 + 4(x)] • N R (<5 C ) - «S 6 (x) 



dN R (6 c ) 
d5 c 



+ 0$(x)], 



(63) 
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where in the second line, we have used a Taylor expansion to first order. The long- wavelength 
"bias at birth" 6*(P) can be defined by the ratio of the correlation functions at the time when the 
objects of scale R are identified, 

bl{R) - (64) 

where we have defined 

i + faW= <*«y» , (65) 

and £ is the usual mass correlation function. The bias at birth is thus given by 

1 dN 

The value of Nr is given from Press-Schechter by equations (15)— (17). Differentiating with respect 
to the critical threshold S c and rearranging gives the following general form for 6*(P): 



6„(fl) = l + i 



P' . 

u(R,t) - 1 



P 



(67) 



where P' = Examples of the ratio — P'/P are given in equations (22)-(22). In the cases we 
consider, the bias at birth will be given by 



K(R) = 1 + — x { 

Or. 



v 2 - 1, Gaussian PDF 

1.45i/ - 1, Texture PDF . (68) 

0.67^ - 1, Non-Gaussian ICDM PDF 



The Gaussian case has been derived previously by Mo and White (1996) and others. The non- 
Gaussian cases are similar to the derivation for v 3> 1 by Kaiser (1984). It is important to note 
that for small v, an exponential PDF such as that from textures or from the ICDM model gives a 
(slightly) larger bias than a Gaussian PDF; but for extremely rare events with v S> 1, a Gaussian 
PDF leads to a substantially larger bias than exponential PDFs. This was previously noted by 
Gooding et al. (1992). 

A similar calculation yields the bias of objects of scale > R to be 

f™dRN R bJR) , 

K > R) = ^ rooj^ ' (69) 

J R dRN R 

As noted by Mo and White (1996), the total bias is simply the average of the individual biases, 
weighted by the mean number density of objects at each scale. An equivalent form for &*(> R) was 
derived by Robinson, Gawiser, and Silk (1998). 

Below, we will primarily consider this bias at birth in the cosmological models we have analyzed. 
This is equivalent to considering models which are dominated by constant merging at the epochs 
we consider, so that objects do not "survive" long enough for their clustering bias to dynamically 
decrease. 
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4.2. Physical Bias from Statistical Bias: Objects at the Jeans Scale 

As far as we know, a necessary condition for matter to be luminous is that the baryons collapse 
to high density. Thus the first criterion we use to explore physical bias is that objects have a mass 
greater or equal to the Jeans mass. 

The Jeans birth bias 6*, j is dependent on the PDF and the degree on nonlinearity at the Jeans 
scale, vj. In Figure 11, we show the value of uj in the ionized and neutral phases, as well as the 
value of vj averaged between the two phases. In most cases, the results are similar to those found 
by OG96: during reheating, the average Jeans scale tracks the nonlinear scale defined by a constant 
value of vj. The early evolution of vj is determined primarily by the power spectrum, and is only 
weakly dependent on the PDF. This is, of course, related to the fact that the Jeans scale depends 
primarily on the thermal evolution of the universe, which we found in the previous section to be 
determined mainly by the power spectrum on small scales. At epochs z < 30 all the models look 
quite similar. This may be because of feedback between star formation and the temperature of 
the universe once reheating has begun, keeping the Jeans scale at a constant nonlinear scale of 
v = 2 ~ 3.5. 

The resulting bias is shown in Figure 12. At high redshift, well before reionization, the Gaussian 
models exhibit the strongest bias, with 6* > 10 at redshifts z > 30 ~ 40. After a pause during 
reheating, during which b* remains in the range of 5 ~ 7, the bias drops quickly after full reionization 
to a level of 6* = 1.5 ~ 2 at z = 3. The T+ACDM model has an intermediate bias, falling to 
6* ~ 10 at redshift z ~ 80, pausing during reheating at 6* ~ 3, and then falling to 6* = 1.5 ~ 2 at 
z = 3. The pure texture T+CDM and ICDM models have a low bias which always remains b* < 3, 
with the T+CDM bias declining slowly to from 2.5 to 1.5 after reionization and the ICDM model 
bias always remaining ~ 1. 

In contrast to many of the observational constraints considered in the previous section, the 
evolution of b* t j depends on both the PDF and the fluctuation spectrum. While the value of vj 
depends primarily on the power spectrum, the relationship between vj and b* t j is quite different 
for different PDFs. 

4.3. Cooling Efficiency and the Origin of Galactic Spheroids 

The Jeans mass is the smallest possible baryonic mass which can collapse, and the value of 
vj derived here gives a lower bound on the average bias of all luminous objects. But a second 
necessary condition which baryons must satisfy in order to become luminous is that they can cool, 
so as to collapse beyond virial density. 

From the cooling grid calculation described in § C and the prescription for the cooling fraction 
e cool equation (C13), we can determine the constraints from cooling efficiency on the formation of 
luminous objects. In particular, at a given redshift z and a minimum cooling fraction e m j n , we 
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can determine the maximum scale R c (z, e m i n ) for which halos collapsing at z have e coo i(i? c ) > e m i n . 
This gives the scale at which at least a fraction e m i n of the baryons in the halo can cool. Given the 
power spectrum this maximum scale can be translated into a value for the degree of nonlinearity 
v c (z, emin)- To summarize, we have the correspondences: 



In Figures 13-14 we show for z < 30 the values of vj derived in the previous subsection along with 
the values of v c (z, e m i n ) for minimum cooling fractions. Each panel shows vj for the same value 
of S7o, with different symbols for the different power spectra and PDF. The lines show contours of 
£min = 10°'~ a5 ' _1 , where all, ~ 1/3, and 1/10 of the baryons can cool. 

What is particularly striking about the results of this calculation is that for z < 10, the ratio 
of the dynamical time to the cooling time for Jeans mass objects is on the order of unity. Thus 
objects at or above the Jeans mass limit are barely able to cool, or can only cool at their centers. 
The implications for galaxy formation are quite profound: because the average cooling times and 
dynamical times are comparable, Jeans mass or larger gas clouds at redshifts z < 10 will on the 
average only have a marginal star formation efficiency The p 1 / 2 ~ r _1 dependence of the ratio 
of the dynamical time to the cooling time means that the centers of these gas clouds will have 
significantly more efficient star formation. 

The next question is whether the hierarchical nature of halo formation leads to gas clouds whose 
luminosity is dominated by orbiting stellar condensations or by a single central stellar condensation. 
In the first case, the gas cooling efficiency of the more massive halos is small enough so that the 
cooling mass M c = M\,e coo \{z , M&) begins to decrease when Mf, is greater than some critical value. 
In the second case, the cooling mass M c must be monotonically increasing with 

We find the second case to be more plausible because the dynamical to cooling ratio is a 
relatively weak function of mass. This fact is straight- forward to understand. At redshifts 20 > 
z > 3, Compton cooling is the dominant coolant at scales of ~ l/i~ 1 Mpc. Because the Compton 
cooling coefficient oc (1+z) 4 , bremsstrahlung cooling does not become dominant until recent epochs. 
Thus at high temperatures, the cooling function A c oc T, and the cooling time t c oc T/A c approaches 
a constant. Since the dynamical time is determined by the redshift, the ratio of the dynamical time 
to the cooling time at high temperature is constant for all masses in a given epoch. The cooling 
fraction thus approaches a constant, so M c oc M and monotonically increases with the total mass 
M. Even with brehmsstrahlung cooling, which rises as A c oc T 1 / 2 , the cooling time would only 
decrease as t c oc T -1 / 2 . The virial temperature is related to the comoving scale by T oc R 2 , the 
cooling fraction will fall as R~ x . Since the mass M increases as M oc R 3 , the total cooling mass 
will still increase as M c oc R 2 oc M 2 / 3 . 

Therefore, the star forming regions of halos formed in the redshift range 10 < z < 3 should 
always be dominated by their central star forming core rather than orbiting stellar condensations 



R c {z,e min ) = max[R], such that e coo i(R,z) > e mi 



(70) 
(71) 




a[R c (z,e min )]' 
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leftover from mergers. Thus a picture of galactic spheroids forming at the centers of much more 
massive gas clouds and dark halos is in place. 

5. Summary and Conclusions 

We have developed a semi-analytic model for the reheating and reionization of the universe. 
Our primary assumption is that UV photons from early generations of star formation provide the 
energy to heat and ionize the primeval gas. We have explored the thermal and ionization evolution 
of the universe in a variety of cosmological models. 

Interestingly, we have found that many of the results found in a set of specific hydrodynamic 
simulations of the A-CDM model (G097) appear to be generic: (1) reheating precedes reionization, 
and clumping is important for understanding both processes; (2) full reionization is a sudden process 
of percolation and occurs near redshift z = 12 ± 5; (3) during reheating, the Jeans mass closely 
tracks the nonlinear mass scale; and (4) hierarchical clustering models are broadly consistent with 
the measured Gunn-Peterson opacity of the IGM. However, models do have significant differences 
between which observations in the near future might be able to distinguish. 

Our major conclusions are as follows: 

1. The most important determinant of the thermal history of the universe for z < 100 is the 
power spectrum at the Jeans scale. The shape of the PDF is only of secondary importance 
to the overall thermal history. For both Gaussian and texture PDFs, models with adiabatic 
CDM-type power spectra, which have a logarithmically diverging fluctuation spectrum o~(R) 
at small scales, lead to reheating beginning at z ~ 30. These we refer to as "late reheating 
models." For models with a power law divergence at small scales, such as textures or ICDM, 
reheating has already begun at z ~ 100. These we refer to as "early reheating models." 

2. In all models, full reionization occurrs as z ~ 10. In models with late reheating, reionization 
proceeds after reheating has plateaued. In the early reheating models, the increased clumping 
of gas caused by earlier structure formation delays full reionization until z ~ 10. 

3. Two important observational constraints which may be able to distinguish between different 
models are the CMBR and the composition of the IGM. With regards to the CMBR, the 
polarization anisotropy caused by Compton scattering off free electrons will be the strongest 
test of when the epoch of reionization occurred. The late and early reheating models give 
different Compton optical depths to recombination which may be testable with experiments 
with sensitivity of ~ 1/iK. With regards to the IGM, the strongest constraint is the baryon 
fraction in stars/quasars. The early reheating models have > 20% of the baryons sequestered 
in stars/quasars by z = 3, while the observational evidence suggests that this fraction is 
smaller. 
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4. We found a correlation between the thermal and ionization history of the universe and the 
formation of luminous matter. In particular, we found that during reheating, the Jeans mass 
scale tracks the nonlinear mass scale in all models, with 5Mj/Mj = 0.4 ~ 1, depending on 
the particular model. After full reionization occurs, the Jeans scale remains roughly constant, 
so 5Mj/Mj grows. 

5. Since 5M/M can be related to galaxy bias, we find that bias at the Jeans scale remains roughly 
constant during reheating, and then drops after reionization. Gaussian models generically 
give much higher biases than the texture and non-Gaussian ICDM models, and is a potential 
discriminant between models. 

We also considered the impact of the thermal and ionization history of the universe on galaxy 
formation. It turns out that Jeans mass objects at z < 10 are, on average, barely able to cool in 
a Hubble time. Therefore, galaxy-sized objects will only cool efficiently in their centers. Galactic 
spheroids, then, may form at the centers of much larger, slowly cooling gas clouds. A rough picture 
of galaxy formation has thus emerged. 

WAC and JPO acknowledge David Spergel and David Wilkinson for helpful comments, and 
Nick Gnedin for providing results from his simulations. WAC acknowledges David Spergel, David 
Wilkinson, and Peter Meyers, who were all on his dissertation defense committtee. 



A. Thermal History of Neutral and Ionized Phases 

The temperature of the neutral (cold) phase will be determined primarily by adiabatic ex- 
pansion and Compton scattering of electrons off the CMBR. We take the electron fraction of the 
neutral phase xj to be equal to the residual electron fraction xq from Peebles (1993) of 

1.2 x 10^ 5 ^ /2 

XI = X0 = Ml b (A1) 

where is the present day value of the total mass density parameter, h is the Hubble parameter, 
and Qt„ is the present day baryon density parameter. We define the gas internal energy and particle 
number density 

3 

u = -nk b T (A2) 
n S 5 5 * (A3 ) 

where kt, is Boltzmann's constant, T is the temperature, n is the total particle density, n& is the 
baryon density, is the particles per baryon, and X is the hydrogen mass fraction. We have 
assumed an adiabatic index 7 a = 5/3. 
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For simplicity, we consider only cooling from hydrogen and electrons, using the compilation of 
Anninos et al. (1998). In erg cm -3 s" 1 , these rates, followed by their reference and formulae (if not 
too complicated) are: 

Compton cooling (heating): Peebles (1993) 

T cc n e = 5.65 x l(T 36 a- 4 (T - n e 
Collisional excitation cooling: Black (1981); Cen (1992) 

r ec n e = 7.50 x 1(T 19 (1 + V 7 ^)- 1 exp (-118348/T) n e n H o 
Collisional ionization cooling: Shapiro and Kang (1987); Cen (1992) 

Ti c n e = 2.18 x 10 -11 k c in e n H o 
Recombination cooling (Case A): Black (1981); Spitzer (1978) 

T rc n e = 8.70 x 10~ 27 T 1 ' 2 r 3 ~ a2 (1 +T6 - 7 )- 1 n e n H + 
Molecular hydrogen cooling: Formula from Lepp and Shull (1983) 

^H 2 cn e 

Bremsstrahlung cooling: Black (1981); Spitzer and Hart (1979) 

F bc n e = 1.43 x 10~ 27 T 1 / 2 [1.1 + 0.34 exp (-(5.5 - log 10 T) 2 /3)] n e n H + 

Here k c i is the H° + e~ collisional ionization rate coefficient taken from the Janev et al. (1987) 
compilation, and T n is the temperature in units of 10 n Kelvin. 

We define the total cooling rate to be 

r to t, c = c (r ec + v ic + r rc + r H2C + v bc ) + r cc , (A4) 

where the clumping factor C is required for the particle-particle cooling rates, but not the Compton 
cooling rate. The equations for thermal evolution in the H I and H II regions are thus 

ui = -5^ur - r to t,c£/"-H (A5) 

u\j = -5-Uij + T ph (l - x u )n H - T tot c x n n H , (A6) 
a 

where the photoheating rate T p h is given by equation (14), and we have assumed that the electron 
density n e is equal to the ionized hydrogen density n H +. The equations for the temperature 
evolution are thus 

hTi = -2-k b Ti - ^-T to t,cXi (A7) 
a 3fi 

k b f n = -2-k b T n + — [r ph (l-x II )-r tot:C x n ], (A8) 

Ob olX 

where X is the hydrogen mass fraction. 

These temperatures are used to determine the comoving Jeans length in each region: 



Rj " m ~ ~ V l2Gp, mH a^ (A9) 
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where kj is the comoving Jeans wavenumber and p is the total mean density (including dark 
matter). The Jeans mass in baryons is thus 

M Ji[m = -jPo,bB?j, [m (A10) 
where po,& is the current mean density in baryons. 



B. Molecular Hydrogen Formation 

We follow the procedure described in Anninos et al. (1998) to solve for the evolution of molecu- 
lar hydrogen. The reactions involving molecular hydrogen are given in Table 2. The formation of H2 
is primarily through the intermediaries H~ and H^, which Anninos et al. (1998) note to be nearly 
in their (small) equilibrium abundances at all times. Neglecting the reaction + H~ — > H2 + H°, 
which is only a small order correction, the equilibrium abundance of H~, x H - = n H -/n#, can be 
written independently of H 2 : 

= k 7 n w x e 

fc 8 n H o + k u n e + /ci 5 n H o + (ki 6 + ki 7 )n H+ + k 23 ' 

where the variables ki are the rate coefhcients with subscripts referring to the reaction number in 
Table 2. Then given n H -, the equilibrium abundance of Hj, x h + = n H +/nn can be written with 
no additional assumptions as 

x + _ kgn u ox R+ + fciira H2 x H + + ki 7 n H -x H + + k 2 AX}j 2 
H 2 fci n H o + ki 8 n e + ki 9 n u - + k 2 5 + k 2 & 

The evolution equation for H 2 is thus 

%h 2 = C — Dxu 2 j (B3) 
where the creation and destruction rates are 

C = kgn u -x u o + kion u +x R o + kign H +x u - (B4) 
D = knn H + + ki 2 n e + ki 3 n H o + k 24 + k 27 + k 2S . (B5) 



C. Gas Cooling in Virialized Halos 

In Press-Schechter theory, the dynamical time for halos collapsing at z co \ is the same regardless 
of its mass, since the virialized density of the halo depends only on the collapse time. Specifically, 
the dynamical time we use is: 



* dyn(Zcol) " V 32Gp 3 J(. col ) (C1) 

Pvir(^col) = A vir [ft(z co i)]p c (z co i), (C2) 
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where p c is the critical density, and the overdensity factor A v i r depends only on Q at the time of 
collapse. 

It is a much more involved process, however, to determine the cooling time i coo i- Although the 
mass density is the same for all halos collapsing at a given epoch, different mass halos have different 
virial temperatures, as well as different temperature histories. Because the primordial chemistry of 
the expanding universe not in equilibrium, we must follow the thermal and ionization history of a 
halo in order to determine its cooling time at virialization. 

Some previous semi-analytic models (e.g., Cole, Fisher, & Weinberg 1994) have simply assumed 
a collisional equilibrium cooling function A C (T) = r to t,c/^H which only depended on temperature. 
Tegmark et al. (1997) considered the minimum mass which could cool within a Hubble time. Their 
results, as well as those of OG96 and G097, showed that molecular hydrogen plays a major role 
in the cooling of early halos (as envisaged by Kashlinsky and Rees 1983 and Couchman and Rees 
1986). The formation of H2, however, is both non-equilibrium and density dependent. Thus, we 
calculate a cooling rate for each halo, assuming a uniform spherical collapse of gas which ends up 
at the virial density. Our cooling function depends on the temperature of the halo and the collapse 
redshift. We do make the simplifying assumption that the same cooling function A c (T,z) applies 
for the entire halo. 

To calculate A c (T,z), we extensively modified the non-equilibrium "TOD" ionization code of 
Abel (1998) so that it modeled a spherical collapse. This code incorporates all of the physics of 
Anninos et al. (1998) and Abel et al. (1998). We assume a uniform sphere of gas and dark matter 
which collapses to one-half its radius of maximum expansion and which heats the gas to the virial 
temperature. Before turnaround, we assume the standard parametric form for spherical collapse, 
given by: 

r = ,4(1 - cos 0) (C3) 
t = B(8-sm6) (C4) 
A 3 = GMB 2 (C5) 

tta = KB (C6) 

rw = 2A, (C7) 

where i t a is the time at turnaround, r max is the maximum radius of expansion, and M is the total 
mass of the halo. After turnaround, we adopt the following simple functional form for the radius 
as a function of time: 

r = A ( 2 - r> + sL^ ) < C8) 

r = ^ (C9) 



For the velocity dispersion a 2 of the gas, we assume shock heating leads to a rise from its value at 
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turnaround <r 2 a to the virial velocity dispersion u 2 ir via the equation 



a 2 = al + (a 2 ir - a 2 J • (l - secli 2 j . (CIO) 

Plots of these functions are shown in Figure 15. 

The calculated value of the ratio idynAcool f° r gas a t the virial density is given in Figures 16 
and 17, for the cases of Qq = 1-0 and f2 = 0.35. Also shown is the comoving Jeans length for gas in 
the neutral phase. The narrow ridge in which the cooling time is longer than the dynamical time 
extending from R rs 0.02, z « 25 to R « 0.01, z « 100 for the f^o = 1 case, and R « 0.04, z ~ 20 
to R ~ 0.02, z ~ 100 for f2o = 0.35, corresponds to when H2 has been destroyed, but atomic line 
cooling is not yet efficient. This "pause" in the cooling efficiency was first identified by Ostriker 
and Gnedin (1996). Note that H2 cooling remains efficient later for = 0.35 than for f^o = 1 (the 
tip of the "banana" shaped region to the left of the ridge extends to lower redshift ). 

This calculation is essentially a simplified version of the Tegmark et al. (1997) calculation, but 
we use the information in the entire "cooling grid" in our work. This is particularly important at 
later times, since the cooling efficiency actually begins to decrease for large masses (this is why we 
have clusters of galaxies rather than gigantic ~ 10 15 M mega-galaxies). These "cooling grids" can 
be precalculated for a given background cosmology, and are assumed to be the same for any PDF. 

Even if the average cooling time of a halo is greater than or equal to the dynamical time, it 
is possible that the central, densest parts of a halo can still cool. To see why, consider a constant 
cooling function A c . Because the cooling time is inversely proportional to the density t coo \ oc p _1 , 
and the dynamical time idyn P~ 1//2 > the ratio of the dynamical to the cooling times is proportional 
to p 1 / 2 : 

^ =^1^1/2 (C11) 

^cool 

If we make the ansatz that only regions cooling on a timescale shorter than the dynamical timescale, 
so that ip > 1, can undergo star formation, then there will be a maximum radius, possibly equal to 
the virial radius, outside of which star formation will not occur. For a singular isothermal sphere, 
the density is proportional to r~ 2 , and the mass within a radius r is proportional to r. Therefore 
the fraction of the baryonic mass of the halo undergoing star formation is simply proportional to 

M 6 [r( y = l)] = r^=l) = /^T (C12) 
M fe , vir r vir V 3p(r) ^ V ' 

If we define ip$ to be the value of (p for gas at the virial density of the halo, then since the mean 
density of the halo p vir equals the local density when r = r v - ir /^/3, the fraction of a halo which is 
cooling faster than the dynamical timescale e coo i(Mfe) will be given by 



e C ool(^4) = min 



V3' 



(C13) 
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This is, of course, an upper limit, since halos will have cores so that the density no longer scales as 
r~ 2 . To account for core radii, we set a e coo i = if (^o/v^ < 0.01, corresponding to a core radius 
equal to about 1/30 of the virial radius. 



D. Calculational Implementation 

D.l. Scaled Equations 

In order to simplify our differential equations (2) and (4), we convert to dimensionless energy 
densities by dividing E and S by eon//: 



£ = 
S = 



E 



S 



We scale the luminosity £(M b ,t,t) equation (27) by eo: 

£(M b ,t,t) 



£(M b ,t,t) 



eo 



yn 



id 



cxp 



We also make the following definitions: 
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x n (l - 



The expression for S in closed form from equation (31) is thus 
S 



n, 



-l 



t - t 



yn / J 



oo ft 

„ , dM b / dtiP(M b , t) ■ N(M b , t, t) ■ C(M b , t, t); 
o Jo 



the expression for £, using equation (43), is 



(Dl) 
(D2) 



(D3) 



(D4) 
(D5) 
(D6) 
(D7) 
(D8) 
(D9) 

(D10) 

(Dll) 



(D12) 



o \ 1/3 ,-oo ,4 

V 4^S ) n ° Jo dMb J di ^ Mb ^ • N ( M *» *' *) • cA/3 ( M b, I t). (D13) 
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Note that although the integrands which determine S and C depend on fu through the selection 
function ip, this dependence is only for /// evaluated at times t < t since C(Mb,t = t,t) = 0. 

These integrals can be simplified by making the definition 
Q(t) = J dR^(R, t) N {oim (R,t)x 

[l-fM^(Rj)^^^-), (D14) 

which is the cumulative source function for star formation, with energy in Rydbergs. The scaled 
source function S is given by smoothing g over a dynamical time, and accounting for the destruction 
of halos through mergers. The smoothing function )C(t,t), from the definition of £, is given by 

D(i) 



JC(t,t) = -rJ— \- — - ) exp 

^dyn V ^dyn 



t - t 



(D15) 



so the scaled source by 



UynJl D(t) 

S(t) = J dtJC(t,t)g(t). (D16) 



Similarly, defining 

g 4/3 (t) = J dRip(R, t)N{ OTm (R, t) 



x 

{ [1 - /,(*)] ««,(* (^1) J*'" , (D17) 

we obtain from equation (D13) 

< = -{-^s) J dt1C(t,t)Q 4/3 (t). (D18) 

In terms of the scaled source function S, the differential equation for /*, the fraction of baryons 
in stars/quasars, is: 

e e ftm H c z 

obtained by differentiating equation (30). The ionization fraction and UV energy evolution equa- 
tions (2) and (4) thus become 

xii = x„[{l-x n )K ci -x n T] + {l-x u )^- (D20) 

J II 

£ = S-[ 1 H+(l-x n )^}£- V - 1 \, (D21) 
having replaced /// by (1 — fn)X. Our constraint equation, derived from equation (43), is 

a=( £n3 

and S and C are defined above. 



A ( - ) , (D22) 
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D.2. Finite Difference Equations 

The differential equations (D20) and (D21) are non-linear and "stiff," so forward differencing 
such as by Runge-Kutta is not stable. Two methods often used in such stiff systems are "implicit" 
and "semi-implicit" differencing (Press et al. 1992). Implicit differencing evaluates the right-hand 
side derivatives at the new location, but is only guaranteed to have a stable solution in closed form 
for linear systems. In our case, the nonlinearity of the equations precludes using purely implicit 
differencing. Semi-implicit differencing linearizes the right-hand side derivatives, inverting a matrix 
at each step, but is not guaranteed to be stable. Furthermore, because of the complicated nature 
of the right-hand side derivatives, the Jacobian matrix required for these methods is extremely 
cumbersome to evaluate. Below, we use a simple "mostly backward" differencing scheme. 

At each timestep, we first evaluate S and ( using equations (D12) and (D13), and update the 
temperature of each phase using equations (A7) and (A8). At very high redshift, S and ( will 
be zero; we leave xjj at the post-recombination residual ionization, and £ and A as zero. When 
S and £ become nonzero, we solve for the new values of A and £ simultaneously from equations 
(D21) and (D22), using the previous timestep value for xjj. We then solve the ionization equation 
(D20) to update the value of xn- Finally, we update the abundance of H2. The finite differencing 
implementation of this procedure follows below. 

In differencing equation (D21), we take the values of (1 — xn) and 77 at the current timestep 
tj. The other values are evaluated at the advanced time step U+i. The finite difference equation 
for dt = tj+i — U is thus 

£ i+ i -£i = dtS i+1 - dt[yH i+1 + (1 - x n )i£i+i] £i+i - ^(Aj+i - A*). (D23) 

Replacing Aj+i with {£ j + i/£j + i) 3 , equation (D22), and collecting terms gives the equation 

(^)\- 1 +£i + iU-V = (D24) 

with 

U = l + dt-yH i+1 + dt(l-x n )iSi+i (D25) 

V = XiV' 1 +£i + S i+1 dt. (D26) 

While this cubic equation has an algebraic solution, it is faster to solve for £j + i iteratively using 
Newton's method. Equation (D24) is of the form 

f(y) = Cy 3 + Uy-V = 0, (D27) 

where y = £i+\. Newton's method iterates to a solution using 

f(vi) Cy? + Uyi - V 
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We use the subscript j so as to not confuse this iterative root finding with our time step labels i. 
The key to Newton's method is a good "first guess" for y. Since U, V, C > 0, we know that the 
solution y satisfies y < (V/C) 1//3 and y < V/U. So our first guess yo = min[(VyC') 1 / 3 , V/U]. In 
cases where U <C C 1 / 3 ^ 2 / 3 or U S> C 1 / 3 ^ 2 / 3 , the iteration will converge in only one or two steps. 
Once a solution for is found, we then use equation (D22) to find Aj+i. 

We use the updated values for £ and A in the ionization equation (D20). We evaluate all values 
at the advanced time step. The finite difference equation is thus 



xn,i+i — x n,i 



dt 



Zff,i+l(l _ XlI,i+l)K c i,i+l - x 2 II i+1 Y i+ i 



+(1 - xn,i+i)Pi+i 



fn. 



(D29) 



Collecting terms gives 

x n,i+idt(F i+ i + K ciii+ i) + xn,i+i 



l + dt[ i+1 Kd ti+1 



fn, 



x n ,i + dt (5 i+ i 



fll,i+i ' 



or equivalently 



with 



A 
B 



x 2 n, i+ i + A x IIti+ i - B = 



1 + dt (f3j + i£j + i/ fn,i+i — K ci)i+ i} 

dt (r m + Ka ii+ i) 
Xi + dt /3i+i£i+ifn,i+i 



(D30) 
(D31) 

(D32) 
(D33) 



dt (r i+ i + Kcij+i) 

This can be solved for xn,i+i either algebraically or iteratively. 

The new values of the electron fraction are used to update the molecular hydrogen abundance. 
Converting the differential equation for the evolution of H2 equation (B3) to a backward difference 
equation gives 

xn 2 ,i + C i+ i dt rr i o A \ 
XH - +1 = l + A+idt (D34) 
where the creation C and destruction D rates are given by equations (B4) and (B5. We use the 
equilibrium abundances of H~ and H^, equations (B1)-(B2). 
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Fig. 1. — Results of reionization calculation with G097 parameters: /// is the filling factor; ( J21) 
is the global average of the ionizing intensity J21; (1 — x) is the neutral fraction (thick: global 
average, thin: ionized region); and (T) is the temperature. The results from the G097 simulation 
are denoted by the stars. 

Fig. 2. — The fraction of baryons in stars/quasars /*, and the clumping factor C for reionization 
calculation with G097 Parameters. The results from the G097 simulation are denoted by the stars. 

Fig. 3. — Jeans and nonlinear mass scales for G097 comparison: solid line marks the volume- 
averaged Jeans mass in baryons derived from our semi-analytic model; the dotted line denotes the 
mass scale at which 8M/M ~ 0.4 at each epoch. The Jeans mass from the G097 simulation, as 
described in Ostriker and Gnedin (1996), is denoted by the stars. 

Fig. 4. — Results of reionization calculation for flat Gaussian models G+ACDM: /// is the filling 
factor; (J21) is the global average of the ionizing intensity J21; (1 — x) is the neutral fraction (thick: 
global average, thin: ionized region); and (T) is the temperature. 

Fig. 5. — Results of reionization calculation for flat texture models T+CDM. 

Fig. 6. — Results of reionization calculation for flat texture models with adiabatic CDM power 
spectrum T+ACDM. 

Fig. 7. — Results of reionization calculation for open Gaussian models. 
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Fig. 8. — Results of reionization calculation for open texture models. 
Fig. 9. — Results of reionization calculation for ICDM model. 

Fig. 10. — Fraction of baryons in stars/quasars for different models. Except for the ICDM model, 
the different line styles denote Q = 1 (solid), 0.45 (dotted), 0.35 (dashed), and 0.25. For the ICDM 
model, which only has one background cosmology, Qq = 0.2. 

Fig. 11. — Degree of nonlinearity v of Jeans scale for different models. The correspondence between 
line styles and are the same as in Figure 10. The upper thin lines are v for the ionized phase, 
and lower thin lines are v for the neutral phase, and the thick lines in between are the average 
values of v. 

Fig. 12. — Bias at birth 6* of Jeans scale objects for different models. The definitions of line styles 
are the same as in Figure 11. 

Fig. 13. — Limits on v from Jeans criterion and cooling efficiency for flat Gaussian and texture 
Models. The adiabatic Gaussian G+ACDM model is denoted by diamonds O, the texture model 
T+CDM is denoted by triangles A, and the texture + adiabatic model T+ACDM is denoted by 
squares □. The contour lines indicate the value of v c (z, e m i n ) with the minimum fractions of the 
mass of the halo can cool in a dynamical time having values of e m i n = 10 ' -0 - 5 ' -1 (solid, dotted, 
dashed lines). 

Fig. 14. — Limits on v from Jeans criterion and cooling efficiency for open Gaussian and texture 
models, and ICDM model. The adiabatic Gaussian model G+ACDM is denoted by diamonds O, 
the texture model T+CDM is denoted by triangles A, and the ICDM model is denoted by the 
crosses x. The contour lines indicate the value of v c (z, e m i n ) with the minimum fractions of the 
mass of the halo can cool in a dynamical time having values of e m i n = 10 0, ~ - 5 ' -1 (solid, dotted, 
dashed lines). 
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Fig. 15. — Top figure gives the radius as a function of time for uniform spherical collapse (thin line) 
and our functional form for collapse and virialization with virial radius equal to half the maximum 
radius (thick line). The lower figure gives the function form of the heating of gas from its velocity 
dispersion at turnaround ofa to the virial equilibrium with velocity dispersion a^ ir . 

Fig. 16. — The ratio of dynamical time to cooling time for Qq = 1, as a function of collapse redshift 
and comoving radius. Also shown is the comoving Jeans length of the neutral phase (white boxes). 
The contours mark where idynAcool = 10 -2 ' -1 '" - 5,0 ' - 5 ' 1,2 . 

Fig. 17. — The ratio of dynamical time to cooling time for 0,q = 0.35, as a function of collapse 
redshift and comoving radius. The notation is the same as in Figure 16 
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Table 1. Cosmological Models Considered and Calculated Results Related to the Thermal and 

Ionization History of the Universe 



Cosmological 






Results 








Parameters 


Reheating 
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early 
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T pp 

ur 
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or late 
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10~ 3 
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= 3) 


Gaussian Adiabatic CDM 














1.0 


0.0 


0.60 


late 


13.1 


10.3 


1.4 


0.069 


0.14 


0.047 


0.45 


0.55 


0.71 


late 


10.6 


7.6 


1.3 


0.065 


0.17 


0.027 


0.35 


0.65 


0.76 


late 


10.2 


7.3 


1.4 


0.068 


0.17 


0.025 


0.25 


0.75 


0.85 
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0.077 


0.19 


0.025 


0.45 


0.0 
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0.17 
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0.35 


0.0 


0.71 
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Textures + CDM 
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5.9 


0.33 


0.11 
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CDM power spectrum 










1.0 
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1.0 
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0.021 


0.45 
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0.47 


late 
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0.35 
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late 


11.6 
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0.016 
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early 
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Table 2. Chemical Reactions Involving H2, , and H 



Reaction number 


Reaction 


(') 


XjO 1 _— , XT— 1 - 

xi -T e — ► Xl "T / 


(8) 


H" + H° -► H 2 + e~ 


(9) 


H° + H+ -> H+ + 7 


(10) 


H+ + H° -► H 2 + H+ 


(11) 


H 2 + H+ H+ + H° 


(12) 


H 2 + e~ -► 2H° + e - 


(13) 


H 2 + H° -> 3H° 


(14) 


H~ + e~ ->■ H° + 2e" 


(15) 


H - + H° -► 2H° + e - 


(16) 


H~ + H+ -► 2H° 


(17) 


H- + H+ -► H+ + e" 


(18) 


H+ + e" -► 2H° 


(19) 


H+ + H" H 2 + H° 


(23) 


H- + 7 -»• H° + e" 


(24) 


H 2 + 7 -► H+ + e - 


(25) 


H+ + 7 -► H° + H+ 


(26) 


H+ + 7 -► 2H+ + e~ 


(27) 


H 2 + 7 -► H* 2H° 


(28) 


H 2 + 7 -► 2H° 
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Fig. 2.- 
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Fig. 10.— 




Fig. 11 — 
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Fig. 12.- 
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Fig. 13.- 
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Fig. 16.— 
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Fig. 17.- 



